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[l] Ten years of CHAMP magnetic field measurements are 
integrated into MFACE, a model of field-aligned currents 
(FACs) using empirical orthogonal functions (EOFs). EOF1 
gives the basic Region- 1 /Region-2 pattern varying mainly 
with the interplanetary magnetic field Bz component. EOF2 
captures separately the cusp current signature and By-related 
variability. Compared to existing models, MFACE yields 
significantly better spatial resolution, reproduces typically 
observed FAC thickness and intensity, improves on the mag- 
netic local time ( MLT) distribution, and gives the seasonal 
dependence of FAC latitudes and the NBZ current signature. 
MFACE further reveals systematic dependences on By, 
including 1) Region- 1 /Region-2 topology modifications 
around noon; 2) imbalance between upward and downward 
maximum current density; 3) MLT location of the Harang 
discontinuity. Furthermore, our procedure allows quantifying 
response times of FACs to solar wind driving at the bow 
shock nose: we obtain 20 minutes and 35-40 minutes lags 
for the FAC density and latitude, respectively. Citation: He, M., 
J. Vogt, H. Liihr, E. Sorbalo, A. Blagau, G. Le, and G. Lu (2012), 
A high-resolution model of field-aligned currents through empirical 
orthogonal functions analysis (MFACE), Geophys. Res. Lett., 39, 
L18105, doi:10.1029/2012GL053168. 

1. Introduction 

[ 2 ] Field-aligned currents (FACs) are important agents for 
energy and momentum transport between the solar wind- 
magnetosphere system and the polar ionosphere-thermosphere 
system. Empirical modeling of the FACs system is required 
by both the ionosphere-magnetosphere coupling community 
and the geomagnetic modeling community [e.g., Lukianova 
and Christiansen, 2006; Venner strom et al., 2007], but is 
challenging due to dynamics and complexity of FACs. The 
basic spatial pattern of FACs was retrieved using TRIAD 
magnetic measurements [Iijima and Potemra, 1976], in 
which FACs follow a double ring pattern with the poleward 
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Region- 1 (R1 for short ) portion and the equatorward Region- 
2 (R2) currents. Through magnetogram inversion techniques 
using ground-based measurements, the FAC system was 
found to be controlled by interplanetary magnetic field (IMF) 
and solar wind (SW) conditions [e.g., Feldstein and Levitin, 
1986], FACs were firstly quantitatively modeled through 
taking the divergence of the ionospheric currents [e.g., Foster 
et al., 1989; Kamide et al., 1981]. The models were con- 
firmed and improved by space-based magnetic measure- 
ments, including datasets of Dynamics Explorer 2 during 
1981-1983 [ Weimer, 2001], Iridium constellation [Waters 
et al., 2001], and Magsat and 0rsted dataset [Papitashvili 
et al., 2002]. For a comparison between different models, 
readers are referred to Stauning et al. [2005]. Taking advan- 
tage of the enormous database and continuous IMF/SW 
coverage, the above models parameterized the global-scale 
FAC geometry and its climatology. However, the models 
suffer from low spatial resolution, characterized by low cur- 
rent density associated with large thickness [e.g., Stauning 
et al., 2005]. This resolution problem arises from the inher- 
ent limitation of the global fitting procedures through which 
the input data are averaged in separate sectors according to 
some ‘governing’ parameters [Anderson and Christiansen, 
2003; Stauning et al., 2005]. Furthermore, many aspects of 
FAC dynamics are not contained in existing models, such as 
the seasonal variation of FAC latitude dependence [e.g., 
Anderson et al., 2008]. In the present work we try to over- 
come the above shortcomings by constructing a Model of 
FACs through Empirical Orthogonal Function analysis 
(MFACE). 

2. Data Analysis and Methodology 

[ 3 ] The CHAMP satellite was launched into orbit with 
inclination of 87.3°, at ^440 km altitude in 2001, decaying to 
~3 10 km at the beginning of2010. Forthis study, ls-averaged 
CHAMP magnetic vector data from ~53,000 orbits collected 
during 2001-2010 are used. The geomagnetic latitude (MLat) 
and geomagnetic local time ( MLT) refer to altitude adjusted 
corrected geomagnetic (AACGM) coordinate. 

[ 4 ] First, we calculate magnetic perturbation vectors 
through subtracting the POMME-6 model [Mans et al, 
2006], both the internal and external parts, from the 
CHAMP data. Then for each auroral oval crossing, we 
identify a reference point of Auroral Current Centre (ACC), 
and estimate a meridian profile of current density by the 
procedures demonstrated in Figure 1 through an example. 
Figure la presents time series of each 6 B component. Each 
of the three components is employed to wavelet analysis 
with a mother function of Paul wavelet of order 3. The 
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Figure 1. Example for the Auroral Current Centre identifica- 
tion and current density estimation, (a) Time series of S B in 
North-East-Center coordinates. 6 B is the magnetic perturba- 
tion vector obtained from CHAMP measurements after the 
subtraction of POMME. (b) The sum (color code) of three 
wavelet spectra for series shown in Figure la, and its integral 
(red line) between 200-1600 km. (c) Estimated current density 
(gray lines), associated with outputs of MFACE (red lines) and 
the model (blue lines) of Papitashvili et al. [2002], In each 
panel, dashed lines indicate the Auroral Current Centre. 


resulting spectra are summarized in Figure lb. The obtained 
spectrum is further integrated in the scale range of 200- 
1600 km, and shown with the red line in Figure lb. The 
maximum in each quarter orbit (0°- max. MLat) is referred as 
the ACC, located at MLat ACC and marked by dashed lines in 
Figure 1. In a 20° -MLat- width interval centered at MLat ACC , 
the perturbations 6 B are used for minimum variance analysis 
(MV A). For the MVA, we define an Orbit-Geomagnetic 
(OGM) coordinate system. In the OGM frame, z is along the 
local POMME magnetic field, x is in the direction of the cross 
product of z and the orbit normal in northern direction, and y 
completes the right-handed orthogonal system. The OGM 
system is an inertial frame for uniform background geomag- 
netic field directions, superior to Geomagnetic coordinates 
defined with geomagnetic East and West [e.g., Wang et al., 
2005] which quickly rotate by ~180° around its z axes every 
time when the satellite goes through maximum MLat. MVA 
in the x-y plane produces the component of 6 B along the 
direction of minimum (maximum) variance. Assuming sheet- 
like FACs, the current density j- is determined as a function 
of MLat according to Ampere’s law, and then mapped to 
1 10 km altitude. Small-scale fluctuations are reduced using a 
Butterworth low-pass filter with a cutoff frequency equiva- 
lent to 220 km wavelength or ~55 km spatial resolution. 
Resulting profiles are interpolated to a MLat grid with an 
origin referring to ACC, termed A MLat grid . 


[ 5 ] In the second step, the gridded profiles are binned in 
each hemisphere for Empirical Orthogonal Function (EOF) 
analysis [cf. Jolliffe, 2002; He et al., 2011], through which 
profiles are decomposed into linear combinations of EOFs, 
j ( y /S.MLat gri d) = j (/S.MLatgrid) + Sj ■ EOFi^AMLatgnd) . 

i 

Here, j is the average of all profiles, s, is the score for EOF t , 
and EOFs are mutually uncorrelated data-derived functions 
arranged in descending order by the variance they account 
for (see also the auxiliary material). 1 The corresponding 
contribution of EOF, to the total variance is ranked in 
Figure 2a. 

[6] After each j : profile is represented with factors of a 
MLat ACC and a group of ,v„ a j z model could be constructed 
through parameterizing the factors. The MLat ACC is regres- 
sed as a function of MLT, Day of Year ( DoY ), IMF clock 
angle 6 IM f in GSM, IMF component in GSM y-z plane B t , 
IMF magnitude B, solar wind speed v S w and AE index. For 
each profile, the parameters refer to the moment when 
CHAMP passes MLat ACC . Mathematically, 

MLat acc(MLT , DoY , 6imf,B,,B, vsw,AE) 

= a 0 + fj ™ a^ or e^ D ° r /365^ 

i= 1 7=1 

3 

+ °‘iMF e 'i 6,m ' + a B,Bi + a$Q(9iMF,B h B.vsiy) 

7=1 

+ cxaeAE + . . . ( interaction and squared terms) ( 1 ) 

Here, i = \J — I , $ is the empirical polar cap potential \Boyle 
et al., 1997], and a’s are coefficients to be determined. We 
exclude outliers characterized by a residual larger than 
expectation at the 5% significance level. Since the IMF 
parameters refer to the Earth bow shock nose, a delay is 
taken into account through replacing 8 IMF , B„ B, v S w in 
equation (1) with 9 IMF {8t), etc. As indicated by the green 
lines in Figure 2b, the best MLat ACC model fit is obtained for 
35-40 minutes lag. When comparing the model predictions 
for local summer with winter, the summer current is located 
at higher latitude around noon by 2-4°, consistent with 
Wang et al. [2005], but at lower latitude around midnight, in 
line with detections of TRIAD and DMSP \Fujii et al., 1981; 
Ohtani et al., 2005]. The green dashed lines in Figure 3 for 
Bz = 5nT indicates an example of this variability. 

[ 7 ] Then, each s t is fitted with the same function as 
equation (1). For the fit, we exclude profiles according to the 
following thresholds: 1) satellite velocity component parallel 
to the MVA current sheet normal less than 3 km/s, 
corresponding to an attack angle larger than ~67°; 2) max- 
imum/minimum eigenvalue ratio (a MVA significance 
measure) less than 2; and 3) angle between the MVA current 
sheet normal and that suggested by the MLat ACC model 
larger than 45°, to ensure the estimated sheet orientation is 
consistent with the statistical R1/R2 rings. Filled portions of 
the boxes in Figure 2a give the fractions of variance 
explained by MFACE. The sum of filled fractions, termed 
Rj, quantifies the quality of total density regression. For the 
first time, the performance of a FAC model is assessed 
quantitatively. Blue lines in Figure 2b show the Rj as a 
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Figure 2. (a) Distribution (the black open boxes) of variance with EOF order, and the corresponding fractions (colorized 
filled portions) captured by MFACE, for Northern (left) and Southern (right) hemisphere, (b) Efficiencies of the FAC density 
model (blue) and of the FAC latitude model (green) as functions of time delay to the 1MF/SW driving at the bow shock nose, 
for the Northern (solid) and Southern (dashed) hemisphere. 


function of time delay to the 1MF/SW. Comparing with the 
R MLat peaks in green, the two Rj peaks in blue at 20 minutes 
suggest that FAC intensity responds quicker than FAC lati- 
tude to solar wind driving. Note that besides the $, we tested 
against other coupling functions for 1MF/SW control such as 
the merging electric field, the Akasofu parameter and the 
Newell function [Newell et al., 2007, and references therein], 
but found $ to be the most effective in maximizing Rj. 
Among geomagnetic indices, the AE index is more effective 
than the Ap, PC, and AL. Since FACs and $ respond non- 
linearly to interplanetary driving [e.g., Anderson and Korth, 
2007, Hairston et al., 2005], MFACE may potentially be 
improved through adjusting equation (1). 

[8] Finally, MFACE is constructed: 

j(MLat, MLT,DoY, 0 IMF ,B h B,v S ir,AE) 

12 

= j(AMLat) + y SjjMLT , DoY, Q lM F,B t ,B, v sw ,AE) 

i= 1 

• EOFi(AMLat) (2) 

Here, A MLat = MLat — MLat ACC ( MLT , DoY, d IMF , B„ B, Vsw, 
AE). As an example, Figure lc presents the modeled current 
density (red), which is much closer to the CHAMP-derived 
profiles (gray) in both intensity and extension than the pre- 
diction (blue) of Papitashvili et al. [2002], For forecast pur- 
poses, we build a model for the AE index using equation (1) 
after replacing AE with B\ h , which is the integral of B s 
(equals to B, when B- < 0, equals to 0 when B- > 0) in the 
preceding hour [A moldy, 1971]. 

3. Results and Discussions 

[ 9 ] Figure 3 presents the polar distribution of FACs at 
moderate conditions as a function of IMF clock angle, sea- 
son and hemisphere, calculated from MFACE. In morphol- 
ogy, all patterns fit the canonical R1/R2 diagram very well 


[Iijima and Potemra, 1976]. In magnitude, the hemispherical 
integrated FAC intensity (top labels) ranges from -4). 7 to 
~3.5 MA in line with existing models [e.g., Papitashvili 
et al., 2002]. The maximum current (bottom labels) ranges 
from ~0.5 to ~2.0 // A/m ~ 2 . more intense than predictions 
of existing models and close to individual events [e.g., Wang 
etal., 2005]. In latitudinal extension, theRl peak spans 2-3°, 
as typically observed [e.g., Iijima and Potemra, 1978]. For 
seasonal dependence, intensities are weaker during local 
winter than summer and equinox, consistent with observa- 
tions from TRIAD and DMSP [Fujii et al., 1981; Ohtani 
et al., 2005], For MLT dependence, the downward and 
upward currents maximize typically in the prenoon and 
postnoon sector respectively, as marked by crosses, in 
agreement with TRIAD observations [e.g., Iijima and 
Potemra, 1978]. For IMF dependence, the most remarkable 
feature is the Ziz-dominatcd size and strength of the FAC 
pattern (also see Figure SI). In particular, under the strongest 
positive Bz (in panels in the top row, marked by the green 
dots as examples), a current pair named NBZ appears pole- 
ward of the R1 sheet, with intensity comparable to Stauning 
[2002], Besides, the By component also influences the FAC 
configuration remarkably. As indicated by the arrows, in the 
Northern hemisphere and under positive (negative) By, the 
R1 sheet on the dawn (dusk) side appears to cross noon and 
extends into R2 sheet on dusk (dawn) side. This By depen- 
dence has been reported before [e.g., Anderson et al., 2008; 
Lu et al., 1995; Stauning et al., 2001; Weimer, 1999], but it 
has not been captured so distinctly in existing statistical FAC 
models. Another .By-related feature that has not been repor- 
ted before is the imbalance between upward and downward 
maximum currents. The dominant peak of the two is high- 
lighted in the plots. In the Northern (Southern) hemisphere 
and under positive (negative) By, the downward current peak 
dominates the upward peak, but is exceeded by the upward 
peak when By changes orientation. 
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Northern Hemisphere Southern Hemisphere 

Jun.Solstice (Summer) Sep.Equinox Dec.Solstice (Winter) Jun.Solstice (Winter) Sep.Equinox Dec.Solstice (Summer) 



Figure 3. Polar distribution of FAC density in the (left) Northern and (right) Southern hemispheres, organized by IMF 
clock angle (marked on left margin) and season (marked on top). AE index for each panel is calculated with the AE model 
mentioned in Section 2. In each panel, numbers in comers show the maximum density (bottom comers, in units of ji Am” 2 ) 
and hemispherical integrated density (top comers, in units of MA), for both the downward (right comers) and the upward 
current (left comers). Maximum upward and downward current peaks are marked by crosses, and are not balanced: the cir- 
cled crosses with shadowed values exceed their counterparts by more than 0.3 /iAm “ 2 in intensity and 20% in ratio. Arrows 
in two colors indicate two kinds of noon time FAC topology. The orientation of By clearly controls the competition between 
upward and downward maximum intensities and the noontime topology, which is opposite in the two hemispheres. In panels 
for Bz = —5 nT, green dashed lines present an example for the comparison of FAC latitudes. In panels for Bz= 5 nT, paired 
green dots indicate examples of NBZ signatures [cf. Stauning, 2002; Vennerstrom et al., 2002], Note that the signature here 
does not cover all NBZ currents because MFACE does not cover the whole polar cap. 


[ 10 ] In Figure 4 we separate the contributions from EOF1 
and EOF2. Overall, EOF 1 plots are controlled by Bz, similar to 
the corresponding panels in Figure 3. In detail, the polarity of 
the double-ring pattern switches sharply around midnight and 
noon as indicated by the black arrows. This differs distinctly 
from the corresponding plots in Figure 3 in which the R1 
currents can stretch into the R2 across noon or midnight. 
Marked by the dotted lines, the midnight switch is identified as 
a signature of the Harang discontinuity, which appears earlier 
(later) under positive than under negative By in the Northern 
(Southern) hemisphere, consistent with previous reports [e.g., 
Rodger et al., 1984]. Different from EOF1, in EOF2 the FAC 
distribution and its polarity are controlled mainly by By. In the 
Northern (Southern) hemisphere, the positive (negative) By is 
accompanied by a triple-sheet structure comprising a dominant 
upward current peak flanked by small downward peaks in the 
morning sector, which shifts to the afternoon sector and 
changes its polarity as By changes orientation. The triple-sheet 
structure in EOF2 could be read as a signature of cusp (also 
termed as Region-0) currents. As reviewed by Cowley [2000], 


the By component influences near-subsolar reconnection, and 
determines the tension and flow of newly-opened field lines. 
For positive (negative) By and in the Northern hemisphere, the 
tension pulls the newly-opened field lines towards dawn 
(dusk), exciting pairs of East-West oriented current sheets in 
the morning (afternoon) sector, and simultaneously in the 
Southern hemisphere pulls field lines toward dusk (dawn) with 
current sheets occurring in the afternoon (morning) sector. The 
current system reaches down into the ionosphere and attaches 
to the R1 ring, thus creating the triple-sheet FAC structure as 
illustrated in Figures 5 and 6 of Cowley [2000]. The dawn- 
dusk asymmetry of cusp current may contribute significantly 
to the imbalance between maximum downward and upward 
current presented in Figure 3. The signature of cusp current is 
stronger during summer than winter, in accordance with pre- 
vious reports [Fujii et al., 1981; Wang et al., 2005], 

[n] In summary, we developed a FAC model MFACE 
from CHAMP data with a novel technique based on EOF 
decomposition. Predictor variables include MLat, MLT, DoY, 
9 /mf, B t , B, Vsw an d AE index. In comparison with existing 
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Figure 4. The same plots as in Figure 3 but for separate contributions from (top) EOF1 and (bottom) EOF2. Dotted lines 
show the By-modulated azimuthal location of the Harang discontinuity. Note that the green lines appear later than the 
corresponding magenta ones in each hemisphere. Black arrows, as examples, indicate that the polarity of the R1/R2 system 
switches sharply in comparison with the corresponding panels in Figure 3. Colored arrows indicate the By-controlled signa- 
ture of cusp currents, where the two hemispheres show opposite polarities. 
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models, the most pronounced advantage of MFACE is the 
enhanced spatial resolution, characterized by more realistic 
FAC intensity, latitudinal extent and MLT distribution. The 
model also captures the season -MLT dependence of FAC 
location. The daytime FACs appear at higher latitude during 
local summer than winter, and the seasonal dependence 
reverses around midnight. Systematic dependences on By are 
manifested in MFACE. Under positive (negative) By and in 
the Northern hemisphere, the dawn (dusk) side R1 current 
extends across noon into the dusk (dawn) R2 current, and the 
upward (downward) maximum current exceeds the down- 
ward (upward) significantly. By modulates also the Harang 
discontinuity which appears earlier for positive By than for 
negative in the Northern hemisphere. All By dependences 
reverse in the Southern hemisphere. Our EOF decomposition 
reveals that EOF 1 represents mainly the /fr-con trolled large- 
scale R1/R2 current pattern whereas EOF2 represents the 
/tv-controlled cusp current signature. In addition, our proce- 
dure allows a rigorous quantitative model assessment and a 
time delay study of FACs in responding to 1MF/SW driving 
at the subsolar bow shock. The correlation analysis indicates 
that the response time is about 20 minutes for FAC density 
and 35^10 minutes for FAC latitude dependence. In brief, 
MFACE captures most of the known features of FAC cli- 
matology, which would contribute to ionosphere-magneto- 
sphere coupling studies and geomagnetic modeling. 
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